Cause and Cure - Deterioration in Accuracy of CFD Simulations 
with Use of High-Aspect-Ratio Triangular/Tetrahedral Grids 


Sin-Chung Chang 
NASA Glenn Research Center, Cleveland, OH 


Chau-Lyan Chang 
NASA Langley Research Center, Hampton, VA 


Balaji Venkatachari 


National Institute of Aerospace, Hampton, VA 


25° AIAA Computational Fluid Dynamics Conference 
Denver, CO, June 5-9, 2017 


This research was sponsored by NASA’s Transformational Tools and Technologies (TTT) project of the 
Transformative Aeronautics Concepts Program under the Aeronautics Research Mission Directorate. 


Outline 


Motivation 
Theoretical Develooment 
Numerical Results 


Conclusions 


AIAA CFD Conference, 2017 


Ss Motivation (I) 


¢ High-aspect ratio triangular/tetranedral elements are 
generally avoided in the near-wall region by CFD researchers 


— Reduced accuracy of gradients 
— Causes Numerical instability 


— Prismatic/Quadrilateral elements can replace high-aspect 
ratio tetrahedral/triangular elements in the near-wall region 


+ Mesh generation becomes complicated (compared to generation 
of pure triangles/tetrahedrons) 


¢ To achieve maximal accuracy, efficiency and robustness, 
triangular/tetrahedral elements are mandatory grid building 
blocks for the space-time conservation and solution element 
(CESE) method 
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Sy Motivation (II) 


With requirements of the CESE method in mind, 
Primary Objective: 


Develop a rigorous mathematical framework that 
identifies the reason behind the difficulties in use 
of high-aspect ratio triangular/tetrahedral mesh 
elements and thereby, helo overcome the 
difficulty. 


AIAA CFD Conference, 2017 


oS Key Outcomes 


e Accuracy deterioration of gradient computations involving triangular 
elements is tied to the elements shape factor 


def . , ’ 
y= sin? Q1 + sin? Q2+ sin? QA3 


Q1,Q@2 and ag = internal angles of the element 


e Degree of accuracy deterioration in gradient computation increases 
monotonically as the value of y decreases monotonically from its 
maximal value of 9/4 (attained only by an equilateral triangle) to its 
minimal limit value of 0* (approached only a high obtuse triangle) 


e Independent of its aspect ratio, the shape factor of right-angled 
triangle has a value of 2, very close to the maximal value of 9/4. 
Hence a grid built from high-aspect ratio right-angled triangles is 
much better for accurate computations of gradients than that built 
from high-aspect-rato and highly obtuse triangles 


e Preliminary extensions to 3D 
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Aspect Ratio of Triangle (1) 


Area A(AP, P2P3) > 0 

For each k = 1, 2,3 

a; = Internal angle associated with vertex P, 

l;, = Length of the side facing vertex Pr 

hy, = Length of the altitude originating from vertex P, 


T > Q1,02,a3 >O0 anda; +a2+a3=7 


1 
i> 0, he > 0, and A= Slx-hy > 0, oe ee 


P, 


Aspect ratio = maximum ratio of edge length to the altitude 


associated with It a 


i; = max{;, _ ee >0 
e |t can also be shown that, the largest side is always associated with the 


shortest altitude 


e This definition is also applicable to ge obtuse triangles (usual 
definition of aspect ratio being ratio of lengths of longest and shortest 
side doesn’t work for such triangles) 
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S& Aspect Ratio of Triangle (Il) 


def 2 
] 2 min — Fe 
(i) n> 5 


(ii) 1] = MWmin ite and only if Q, =A. =—=a3 = = 


w | > 


(iii) lim n = +00 


min(a, ,2,a@3)>0TF 


and (iv) 7>>1 if and only if min{a;,a2,a3}<1 


Aspect ratio 7 
e Attains its minimal value when the triangle is equilateral 


e Attains a large value when one of the internal angles of the 
triangle has a very small value 
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S& Shape Factor of a Triangle 


shape factor, y, of a triangle is defined as 
ry eS sin? a, + sin? ag + sin? a3 
It can be shown that: 


9 
22 
Oe 


9 
(ii) ¥=Ymae 7 a1 = 02 = 03 = 5 Hh =b=1>0 Fi 
(iii) y > OT & max{aj, a2,a3} 3 a & min{aj, a2,a3} 3 0T > 7 4 +00 
and (iv) if one of {a1,a2,a3} = 5 =>7y=2 (even if one of {a1,a2,a3} < 1 and thus 7 > 1) 
Note: 
<> means “if and only if” 


=> means “implies that” 


e Shape factor of a high aspect ratio right-angled triangle (y = 2) 
not far from that of an equilateral triangle (y =9//4) 
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Sy Aspect Ratio and Shape Factor of a Triangle 


P, P 
— at a -'N ae 
P, 5 P,P, P, 
0O<7-a; <land0<a2,a3 <1 a, a3 x%7/2and0<ag<1 
(a)n >> landy<1l (b) 7 >> landy2 


e Both triangles shown above have high aspect ratios (7), 
because one of their internal angles is very small ( << 1) 
- But their shape factors (y) are vastly different 
e Because the value of ai is close to z and therefore az, a3 <1, 
the triangle depicted in left (a) has a high aspect ratio, and is 
also highly obtuse and thus has a very small value of shape 
factor 


- Accuracy of gradient computation from this triangle will be 
poor in comparison to the one from right (b) 
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S Gradient Vector and Directional Derivatives 


A(AP, P2P3) = }d| /2 > 0 
where, 


def 


YQ2—2X, Yo2- Yi 
Y3—2%1, Y3- Yi 


= (2 — £1)(y3 — y1) — (%3 — L1)(y2 — yi) #0 


Py (x; ¥) 


For each k = 1, 2,3, let 


(i) d, be a scalar value assigned to the spatial point Px, 
and 


(ii) P, denote a point in the x — y — ¢ space with 
y coordinates (Xx, YR, Pk): 
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S 


Gradient Vector and Directional Derivatives (Il) 


It can be shown that points Py, Pp, and P3 lie on 
the plane in the x — y — @ space defined by 


o=anr+by+ec 


where, 
ZY. Pi 
d2- 91 Yo- Yi 221 O7= O71 LQ Yy2 2 
3-91 Y3-Y1 L3—-X1 3-91 r3 Y¥3 «3 
a= a a = ; , and c= ‘ 


Thus the x— and y— components of Vo on the plane I’ are 
v, = $ =aand y= 3 =b (on [) 
respectively. 
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S Gradient Vector and Directional Derivatives (Ill) 


Given ns two vertices P; and P;, the directional derivative [;,; 


FP 
along Py, Pr, is defined by . ] 
1 
i. ee rene tp P, 
24g 
l3 
where, P, 
—— 
sig © [PLP)| = of (es — 24)? + (yj — ws)? > 0, 1A j and i,j = 1,2,3 

we have 

4 = 52.3 — 63:23 lo = 33.198 1.3 and /3 = 1.9 821 
and 

7 _ o2- 91 b3 — o2 b1 — $3 

—2.1 = H1,2 = 


» 13,2 = 12,3 = and — i383 = Pet — 
ls ly ly 
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Sy Gradient Vector and Directional Derivatives (IV) 


From earlier equations, we have 


l3-pratlh - pa3 +le- p31 = 0 


and 
Ha + 11,2 = 3,2 + 2,3 = 11,3 + U3,1 = 0 


e Six directional derivatives of @ are linked by four 
independent conditions 


e Any one of them can be determined in terms of any two 


independent directional derivatives associated with two 
different sides of AP, PP 
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S 


Gradient Vector and Directional Derivatives (V) 


(i) For each k = 1, 2,3, let Ag, denote the variation (or numerical error) of @, 


at some given time level introduced through a time-marching procedure. 


(ii) Let Av, and Av, denote the variations of vz; and v, (w— and y— 
components of V¢) induced by the presence of Ad;,, k = 1, 2,3. 


Adz — Adi yo-yi rg—-2, Ado—Ad 


Ads — A — = Ada — A 
he 03 - ys Yl g Av = tg 2i ~ P1 


(iii) For any pair of 2 and j with i #7 and i,7 = 1, 2,3, let Ay;,; denote the 
error of f4;,; induced by the presence of Ad; and A@,;, respectively. 


Then, 
Ad; — Agi 


—ALg,5 = Api = 
U5] 


lz - Apia +h - Apa,3 + lo - Ap31 = 0 


and 
Apes + Api = Apz.2 + Ape3 = Apiy3 + Auzi = 0 
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Sy Gradient Vector and Directional Derivatives (VI) 


Let, 
e = Ap23, €2 s Apuz3.1 and €3 s Api,2 
Then, it can be shown that: 
(i) ly -€, tlg-€2 +13 -€3 =0 
and 


(ii) For any pair of 7 and j with i #7 and i,j = 1,2,3, Ap,;,; can be uniquely 
determined in terms of any two of €1, €2,€3 assuming (rz%, yx), k = 1,2,3 are 
given. 


Moreover, with the use of law of sines, 


(sin a1) €, + (sina) €2 + (sinag) €3 = 0 


The above equation links the errors €1, €2 and €3 of the directional derivatives 
evaluated along the three sides of AP, P>P3 
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Error Amplification Factor (1) 


Let 


pit {Ave + Bus gle)? +(e)? + (es)? £0 


[(ex)* + (€2)* + (€3)?]/3 


(i) By definition, R is square root of the ratio of the two simple averages, i.e., 
the simple average of (Av,.)* and (Av,)*, and that of (€,)?, (e2)* and (e€3)? 


(ii) \/[(Av,)? + (Av,)?]/2 = error norm associated with evaluation of Vé 


(iii) \/[(e1)? + (€2)? + (€3)?]/3 = error norm associated with evaluation of 
directional derivatives along the three sides of AP, P2P3 


(iv)As such R is an error amplification factor, measuring how large 


the error norm for evaluating the gradient vector Vo is amplified 
from that for evaluating the directional derivatives along the three 
sides of AP, P>P3 
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Error Amplification Factor (II) 


(i) The value of R is a function of any two chosen independent parameters (say 
€, and €2) among €1, €2 and ¢3. 


(ii) Jo_(7) < R< Vor+(7),0<7< 9/4 


where, 
def 
os(7) = % [3 42,/0/4 —7|,0<7< 9/4 
and (iii) Each of the two bounds \/o_(y) and ,/o1(y) can be attained 
by R with some special ratio between ¢, and €». 


Thus \/o_(y) and ,/oi(y) are the greatest lower bound and the least 
upper bound of Rf, respectively. 


Note: y is the shape factor of AP, P>P3 
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S 


Properties of o_(y) and o1(y) and Their Implications 


As the value of y decreases from its maximal value 9/4 to its minimal limit 07, 


+ 


(ii) the value of o_(y) decreases monotonically from 1 to (5)*. 


1 
2 
As such one can show that : 


(i)o_(V)=loeop(y)=ley= 2 <> AP, P2P3 is equilateral 
(ii) § <o_(y)<1l<o1(7), if0<y7< 4 


(iii) o-(y) > (6)? S 04 (7) 9 +0 & 7 3: 0F & marf{at, a2, a3} > 7 


and 


(iv) max{q1,02,a3} = 4 >y=2> /Joily) = Vix 1.225 
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Sy Properties of o_(y) and o1(y) and Their Implications 


P, P 
p p 
2 3 P, Pp 
ea 3 
0O<7-a; <land0<a.,a3 <1 a, a3 x7/2and0<ag<1 
(a) AP, P2P3 is highly obtuse with (b) AP, P2P3 is not obtuse with 
n>Ilandy<l n>1landyr2 


(i) R=1 for any (€1, €2,€3), i.e., no amplification error for gradient 
computation, if and only if AP, P)P3 is equilateral 


(ii) A very large error amplification for evaluating Vo could occur if 
AP, P2P3 is highly obtuse 


(iii) Only minor error amplification for evaluating Vé would occur if 
AP, P2P3 is nearly a right-angled triangle 
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Sy Extension to 3D (for tetrahedrons) 


¢ Mathematical framework presented here for triangles can be extended in 
a straightforward manner for tetrahedral elements. 


— Algebra becomes complicated 


¢ For a regular tetrahedron (i.e., one with all four of its faces being 
equilateral triangles), it can be shown that R =1 for all possible 
combinations of the numerical errors associated with the directional 
derivatives evaluated along all six edge-directions of the tetrahedron 


— No accuracy deterioration in gradient evaluation for a regular 
tetrahedron 


¢ For a tetrahedron in which three right internal angles share a common 
vertex, it can be shown that the least upper bound of Ris V2, a number 
slightly larger than 1. 


— Mild accuracy deterioration in gradient evaluation 
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Numerical Results from CESE 


S& Propagation of plane acoustic wave 


* Acoustic wave sent through a 


viscous-type mesh with non- 


uniform spacing Near-Wall Stretched Mesh 
¢« Preserving the phase and amplitude 
of the wave without damping ona 
non-uniform mesh is challenging 


Largest Aspect-ratio location: 


bottom of the domain =225 6 07 0D 03 
Density Contour 
The time-accurate local time | i ae 


stepping procedures, along with 0.99992 0.99995 0.99998 1.00001 1.00004 
ability to minimize dissipation in el 
CESE helps with this. Time-accurate local time-stepping 


Constant time step 


No attempt to control 
dissipation/CFL 
variation 
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S& RANS computations of Mach 2 flow over adiabatic flat plate 


5 Mach 2 fiat plate adiabatic wa 


Mach 2 freestream, Re =15x10®, Adiabatic Wall 


Reynolds-Averaged Navier-Stokes (RANS) —— 
computations 0.0085 Citad $A Model 
Spalart-Allmaras (SA) and Mentor’s Shear-stress 9: 
transport (SST) models B a ooe 


Triangular Mesh (105,000 elements) with highest 
aspect ratio of 3000. 


Mesh (the structured mesh was sliced into 
triangles) and comparison data obtained from 0.001 
NASA Langley turbulence modeling resource 0 05 
website 


0.002 


0.0015 
Coefft of skin-friction along plate 


x (m) 


Al Re_theta = 10,000 
e74d SST Model 
Cfi3d, SST-V Model 
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Velocity profile comparison at a location where «——_ * 
Re (based on momentum thickness) =10,000 ° 


10 
y+ = normalized wall distance = y * u_tau / kinematic viscosity 
U + = normalized velocity = u/ u_tau 
u_tau = friction velocity = sqrt( Shear stress at wall/density) 


10 
log(y+) 
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Ss Hypersonic flow over a 15° — 25° double cone 


Mach 11.3 Laminar flow - 
¢ CUBRC experiment (Run 35) 
313,000 Triangular elements, Max aspect ratio: 1400 


Triangular/tetrahedral elements traditionally avoided 
in hypersonic computations (carbuncle, poor heat- * 


transfer prediction) 15 
Separation bubbles around the two corners and heat 
transfer coefft predicted well. as ° 
small discrepancy in heat transfer peak levels being 


investigated upon through grid refinement : ee es 


CESE Solution 


3 
Oo Experimental Data 


CESE Solution °° 0.1 = heat transfer 
2.5 oO Experimental Data 
coeftt 
_, 0.08 fr ? : 
2 . 
Z 0.06 a 
< e? 
aS} 8 
5 ; Py 
= 0.04 ; 
Phy 


v 


ag 
ro) 
iN) 


oO 
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Sy Summary 


¢ Discussion on fundamental accuracy issues when a mesh with large- 
aspect -ratio triangular/tetrahedral elements are used in CFD simulations, 
especially inside the boundary layer. 


* Sources of inaccuracy related to triangular grids was identified 
theoretically 


— closely related to the elements shape factor (and not simply aspect ratio) 


— deterioration in accuracy of gradient evaluation can be mitigated by use of 
right-angle / equilateral triangles for gradient evaluation 


¢ Effectiveness of CESE in handling high-aspect ratio triangular meshes 
demonstrated through a few viscous flow computations 


- Gradient evaluation procedure take advantage of the analysis shown and 
utilize triangles with good shape factors for gradient evaluation 


« Extension of mathematical framework to tetrahedral elements is in 


progress. 
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Backup Slides 


S& Original approach of evaluating derivatives in CESE 


S 


Original approach of evaluating derivatives in CESE 


S& Original approach of evaluating derivatives in CESE 


S& Original approach of evaluating derivatives in CESE 


Triangles with bad shape factor that are used in gradient 
computations 


S& New edge-based derivative approach 


Triangles with good shape factors to be used in gradient 
computations, obtained by adjusting stencil (can also be 
used to adjust dissipation) 


